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Abstract 


The  paper  addresses  results  of  numerical  experimentation  with  an  a-posteriori  error  estimator 
which  was  theoretically  analyzed  in  [1].  The  paper  concentrates  on  the  problems  of  the  selection 
of  benchmarks  which  is  directed  to  a  numerical  verification  of  basic  theoretical  features  of  the 
estimator. 


1  Introduction 


During  recent  years  much  interest  has  been  focused  on  the  design  of  a-posteriori  error  esti¬ 
mates  in  the  finite  element  method,  their  experimental  verification  and  their  use  for  adaptive 
procedures.  We  refer  for  example  to  [1-24]  and  further  citations  there. 

The  a-posteriori  error  estimators  are  often  derived  by  purely  heuristic  reasoning  and  without 
mathematical  analysis.  If  a  mathematical  analysis  is  made  it  adresses  most  often  only  asymp¬ 
totic  properties  and  special  meshes.  Usually  the  estimators  are  numerically  analyzed  (verified) 
on  a  set  of  examples  (benchmarks)  which  are  selected  at  best  by  heuristic  reasoning.  Often 
these  benchmark  computations  are  used  for  the  design  of  a  correction  factor  which  is  used  to 
improve  the  quality  of  the  estimators  for  this  set  of  benchmarks. 

This  pragmatic  appro2K;h  to  test  and  compare  the  estimators  by  benchmark  computations  has 
serious  shortcomings.  The  quality  of  an  estimator  is  usually  very  sensitive  to  the  structure 
of  the  solution  and  the  meshes  used.  The  performance  of  an  estimator  can  be  very  different. 
It  may  depend  on  whether  the  meshes  are  nearly  translation  invariant  or  general,  or  whether 
the  meshes  are  general  or  such  that  the  error  indicators  are  nearly  equal  on  all  elements  (i.e. 
the  meshes  are  equilibrated),  or  whether  the  solution  is  smooth,  has  point  singularities  or  is 
generally  unsmooth,  etc.  Some  of  these  aspects  are  adressed  in  [5].  Hence  the  benchmark 
computations  could  motivate  misleading  conclusions.  Each  benchmark  computation  should  be 
a  representative  of  a  more  or  less  precisely  defined  class  of  computational  problems  and  the 
conclusions  used  for  this  class  of  problems. 

In  this  paper  we  will  address  various  aspects  of  the  relation  between  theoretical  understcuiding 
of  an  estimator  and  the  benchmark  selection.  We  will  adress  a  certdn  estimator  which  is  used 
in  practice,  and  discuss  its  properties.  We  do  not  intend  to  compeure  it  with  any  other  error 
estimator. 


2  The  Model  Problem 


Let  n  €  be  a  bounded  polygonal  domain  and  let  F  be  its  boundary.  Assume  further  that 
17  =  (J]^i  JJi,  where  fl,  is  a  polygonal  domain  (with  boundary  dfli)  and  that  Hi  D  fij  =0  for 
i  ji  j.  We  will  be  interested  in  solving  the  following  problem: 

du 

-div(aVu)=:/  in  D,  u  =  0  on  Fo,  a-^s^gonTfi,  FoUFat^F,  (2.1) 


where  n  is  the  outward  normal  to  fl.  We  assume  that  /  €  9  ^  o  >  oo  >  0 

is  constant  on  flj,!  =  We  will  assume  that  either  meas(-Fp)  ^  0  or,  if  F^)  =  0, 

JfiJdxdy  4-  J^gds  =  0  and  u[A)  =  0,  A  €  IT.  (Because  of  the  assumptions  about  /  and  g, 
u  €  H^{d),a  >  1,  and  so  u(A)  =  0  makes  sense.) 


We  will  consider  the  weak 
space  and  let 


solution  of  (2.1).  To  this  end  we  denote  by  the  usual  Sobolev 

=  {u  e  ^‘(n),«  =  0  on  Fz,} 


Further  we  let 


(2.2) 


1 


(2.3) 


be  the  bilinear  form  defined  on  Hg{Q)  x  //o(^) 

F(v)  =  /  fvdxdy  +  /  gvds. 

J{1  Jv  jv 

be  a  linear  functional  on  Ho{Vl).  Then  the  weak  solution  u  of  (2.1)  in  Hq  exists,  is  unique  and 

B{u,v)  ^  F{v)  'iv  €  (2.4) 

If  To  =  0,  we  consider  H^(Q)  instead  of  /fo(^)  assume  u(/4)  =  0,  which  is  well-defined,  as 
we  have  in  fact  higher  regularity  for  u. 

Let  us  consider  a  family  M  of  triangular  meshes  M  on  Cl  satisfying  a  uniform  minimal  angle 
condition.  We  denote  by  T  the  closed  triangles  of  the  mesh  M,\JT  =  Tl.  We  will  assume  that 
ifT  £  M  then  interior(T)  C  flj  for  some  /,  i.e.  T  lies  in  a  domain  where  a  is  a  constant.  To 
every  mesh  M  we  associate  a  parameter  t{M)  >  0,  for  example  t{M)  =  where  N  is  the 
number  of  the  nodes  that  do  not  belong  to  To  or  t{M)  =  h  where  h  =  max(d(T))  and  d{T) 
denotes  the  diameter  of  7}.  We  will  always  assume  that  h  =  max(d(7^))  — ►  0  as  t{M)  —*  0. 

Further  let  S{M)  be  the  set  of  all  continuous  functions  on  Cl  whose  restriction  on  a  T  £  M 
is  a  linear  function.  Obviously  5(M)  C  H^(Cl).  Let  So{M)  =  Hq{CI)  fl  S{M)  and  ^^(Af)  = 
{u  €  5(A/),  it(/l)  =  0}.  Then  the  finite  element  solution  us  £  So{M)  satisfies 

B{us,v)  =  B{n,v)  Vr  €  5o(M)  (2.5) 

where  u  is  the  exact  solution  of  (2.1).  If  To  =  0,  we  replace  5o(M)  by  Sy4(M).  We  will  consider 
a  family  of  finite  element  solutions  us  associated  to  M  and  assume  that  ||u  —  us(a/)I|£  = 
B(u  —  as(Ar)*w  —  -♦  0  as  t{M)  -*  0.  By  (2.5)  we  associate  to  any  solution  u  its  finite 

element  solution  us  =  us(A/,u)  and  an  error  estimator  €{u,M)  that  will  be  defined  in  the 
next  section.  In  fact  £  depends  on  us  and  M  but  because  of  (2.5)  we  write  £{u,  M).  This 
estimator  approximates  the  energy  norm  of  the  error,  i.e.  £{u,  M)  as  ||e||£;  =  {B{e,  e))'/*,  where 
e  =  u  —  Us. 

3  The  Error  Estimator 

Let  M  £  M,  Ti,Ti  £  M  he  two  triangles  with  common  side  7(.i  as  shown  in  figure  (3.1).  By 
n(,i  we  denote  the  outer  normal  from  Tj  to  Tj;  we  will  also  write  7/  =  7i„  and  7i  =  T^ut,  as  Tm 
is  the  triangle  for  which  the  normal  nj,i  is  the  outward  one. 

Let  us  now  define  for  every  side  7/,,- 

^  («'^«s)  Ixoat  •  "M  -  (oVus)  |7i„  •  nj,.  (3.1) 

the  jump  of  across  the  side  7j,i.  This  value  is  independent  of  the  choice  of  n(,i,  i.e. 

For  each  triangle  T  £  M  let 

=  —  [div(aVu)  It  -I-  div(aVus)  |t]  +  /  It  .  (3.2) 
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Figure  3.1:  Scheme  of  neighbouring  elements 


We  may  write  for  any  v  €  J/q  residual  equation 


+  ,5/  I°^I  "* 

where  G  stands  for  the  set  of  all  interior  edges  of  M  and  n  for  the  outward  normal  to  fl.  (If 
Fo  =:  0  the  equation  (3.3)  is  the  same,  but  v  € 


For  any  T  €  M  yte  define 


^r  =  :^^J^Rrdxdy 


where  |T|  is  the  area  of  T.  (Hr  is  the  L^{T)  projection  of  the  interior  residual  onto  the  set  of 
constants.)  Further  for  each  7  €  Fat,  where  7  is  a  side  of  a  triangle  T,  we  define 

where  I7I  is  the  length  of  7.  (11,  is  the  .^*(7)  projection  of  the  boundary  residual  onto  the  set 
of  constants.) 

For  each  edge  7  of  the  mesh  let 

[.^1  ifiec 

J,  =  <2II,  if7Crjv  (3.6) 

0  if7C  Fd 

and  for  each  triangle  T  €  M  we  define  the  error  indicator  Tfr' 


which  approximates  the  energy  error  on  T.  Further  we  define  the  error  estimator 

/_ 


e{u,M)  :=  f  I]  qr) 

\T€M  / 


where  x  is  a  suitable  constant  which  will  be  defined  later.  It  is  convenient  to  write 

£*(u,M):=f  (3.9) 

\T£M  / 

for  the  choice  a  =  1.  Let  us  note  that  the  error  indicator  is  meaningful  for  more  general  / 
and  g  than  we  have  assumed.  Nevertheless  our  generality  is  sufficient  for  practical  purposes. 
Obviously  the  set  of  problems  we  are  discussing  is  precisely  described. 
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4  The  Quality  of  the  Error  Estimator 


In  Section  3  we  defined  the  error  estimator  £{tt,  M)  which  approximates  |1c|1b.  In  the  following 
we  will  discuss  the  relation  between  f (u,  M)  and  ||e||e,  especially  the  effectivity  index  of  the 
estimator 

.  £(u,M) 

Mb  ■ 

Before  addressing  this  problem  we  introduce  some  notation  which  we  will  need  in  the  following. 
Let  5’’(M)  be  the  space  of  all  continuous  functions  on  0  whose  restriction  on  T  €  M  is  a 
polynomial  of  degree  p.  Let  So(M)  =  5”  D  “s  €  So(M)  be  the  finite  element 

solution  of  problem  (2.1),  i.e. 

5(u5,w)  =  B(u,t;)  Vt>e5*(A/) 


(If  To  =  0,  we  use  5^(M)  instead  of  So(M).)  Obviously  U5  is  the  p-version  solution  of  problem 
(2.1).  For  more  about  the  h  —  p-version  we  refer  to  (6).  We  have  for  any  A:  >  1: 


1I«S  -  lilU  <  Cp-<*-‘>||u||„»A“"<'’-*-» 

(4.1) 

with  C  dependent  on  k,  but  not  on  p  and  h. 

Remark:  Although  in  (4.1)  the  maximal  element  side  is  used,  an  analogous  element  by  element 
estimate  holds. 

Now  we  are  able  to  state  the  basic  theorems  about  S*(u,M)  .  (See  [1]  for 
and  Section  5  for  an  outline  of  the  main  ideas.) 

the  detailed  proof 

Theorem  4.1  For  any  integer  p  there  is  a  constant  depending  only  on 

ar  of  the  triangles  used  and  on  p,  such  that 

the  minimal  angle 

Me  <  K^e^iu,  M)  +  |1«  -  tillls  +  0{hi) 

(4.2) 

and  the  constant  K*  grows  logarithmically  with  p. 

Theorem  4.2  There  is  a  constant  Ki  >  0,  depending  only  on  ar  and  on  the  jumps  of  a,  such 
that 

n«.M)</f,||e||e-|-0(Ai).  (4.3) 

The  constants  K\  and  K*  are  computable.  We  have 

(4.4) 

and 

-  (m^Cj-)<?., 

(4.5) 

where 

Q.  =  max  (-^T^). 
rrrt^  a\r> 

For  2  <  p  <  8  we  have  the  estimates 

(4.6) 

.548logi  p8in"i(^)  <  Cf  <  .8131og^  p8in“i(^) 

(4.7) 
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a 

c;. 

C^T 

p  =  2 

p  =  4 

p=^6 

p  =  8 

7.5 

.051 

2.390 

3.306 

3.660 

3.988 

15.0 

.072 

1.682 

2.341 

2.609 

2.839 

22.5 

.087 

1.363 

1.918 

2.156 

2.343 

30.0 

.099 

1.169 

1.670 

1.895 

2.058 

37.5 

.108 

1.035 

1.508 

1.727 

1.876 

45.0 

.115 

.939 

1.400 

1.615 

1.757 

52.5 

.119 

.875 

1.334 

1.547 

1.684 

60.0 

.121 

.850 

1.309 

1.522 

1.657 

Table  4.1:  Constants  and  C'-j-  for  various  angles  or 


3.45sin-i(^)  <Cr<  5.85sm-i(^)  (4.8) 

where  ar  is  the  minimal  angle  of  T.  In  Table  4.1  we  give  the  values  of  and  Cj. 

Let  us  note  that  for  h  small  only  very  few  elements  have  a  side  on  the  interface  9fl;  we  can  in 
fact  use  Qa  =  1. 

Remark:  We  have  dealt  only  with  the  Laplace  equation.  In  this  case  we  use  the  minimal  angle 
aj  of  the  triangle.  For  a  general  differential  equation  we  have  to  replace  it  by  which  is  the 
minimal  angle  of  the  triangle  after  a  local  transformation  of  the  general  equation  to  Laplace’s 
equation  (this  is  always  possible). 

The  term  0{h^)  in  (4.2)  and  (4.3)  depends  only  on  ||/||//>(n)  ll5llK>(n> 

emy  u  the  term  ||«  —  uj||£  can  be  made  arbitrarily  small  by  (4.1)  using  sufficiently  large  p.  We 

note  that  the  smoothness  of  u  is  governed  by  the  smoothness  of  /  and  g,  transition  of  boundary 

conditions  and  smoothness  of  F.  It  is  unessential  that  the  estimates  (4.2)  and  (4.3)  have  an 

asymptotic  char2w:ter.  For  example  if  /  =  const  on  every  T  and  g  —  const  on  every  triauigle 

side  of  Tff,  the  term  0{h^)  disappears.  Hence  as  t{M)  —*  0,  the  term  0{h^)  is  negligible  in 

general. 

Theorems  4.1  and  4.2  indicate  that  the  error  estimator  F*(u,  M)  can  be  either  smaller  or  bigger 
than  the  true  error.  Further  the  error  estimator  loses  accuracy  when  or  — ♦  0  and  for  unsmooth 
u  the  reliability  of  the  estimator  decreases  because  a  higher  p  is  needed  in  the  estimator. 

The  coefficient  k,  which  scales  the  error  estimator  {€{u,M)  =  k€*(u,M)  ),  can  be  calculated 
in  different  ways.  We  can  use  the  geometrical  average  of  the  bounds  in  (4.4),  (4.5).  Another 
possibility  is  to  choose  k  such  that  the  error  estimator  is  asymptotically  exact  for  uniform 

meshes  of  equilateral  triangles  and  /  =  0.  In  this  case  we  get  k  =  (^8^3)"^  «  0.26864  which 
will  be  used  for  all  examples. 

Theorems  4.1  and  4.2  indicate  the  expected  theoretical  properties  of  the  estimator  well.  The 
structure  of  Theorems  4.1  and  4.2  clearly  suggests  the  set  of  benchmark  experiments  that  should 
be  performed. 


5 


5  The  Outline  of  the  Proof  of  Theorems  4.1  and  4.2 


For  simplicity  assume  that  /  and  g  are  constant  on  each  7  resp  /  €  Fyv.  (The  general  case  is 
disctissed  in  [1].)  With  e^  :=  uj  —  U5  €  Sq,  we  have 


lletlfi  <  ll«i.ll£:  +  llw  -  UsHe  (5.1) 

and 

l|c,|||  =  B(ep,  Cp)  =  B(u5  -  u,  Cp)  +  B(u  -  us,  Cp)  =  B(e,  Cp)  (5.2) 

For  any  continuous  function  v  let  denote  the  line<ir  interpolant  of  v  on  7.  Using  (3.3),  we 
have: 


llcpill:  =  5(e,cp)  = 


< 

< 


E 


ICdT- 
.  2 


nT(ep  -  e[)dxdy  +  ^ 


’It 

TeM 

53  »?tCJ||cp1|£,t 

Tew 


IC9T 


where 


=  sup 


/ ml  /t(w  -  +  5  E/cst  Illy  /;("  - 

/riVH^dxdy  j 


1/2 


Here,  Pp  and  Pa  denote  the  spaces  of  all  polynomials  of  degree  p  or  0  on  T,  respectively. 
Hence 


with 

K*  =  sup  Cf . 

Tew 

Using  (5.1)  we  get  (4.2)  without  the  term  0(h^^^).  This  is  related  to  the  simplification  assump¬ 
tion  about  /  and  g  made  in  this  section. 

Let  us  now  address  inequality  (4.3).  For  any  triangle  7  as  shown  in  Fig.  5.1,  let  cpo  be  the 
cubic  bubble  function  vanishing  on  67  with  ipo(Qo)  =  1,  where  Qo  is  the  barycenter  of  7.  Let 


Figure  5.1:  Notation  for  triangular  elements 
further  T^{7)  be  the  space  of  quadratic  functions  vanishing  at  the  vertices  of  7  and  let 
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be  the  canonical  basis  of  this  space  (i.e.  <Pi{Qj)  =  =  1,2, 3, j  =  1,...,6)  and  let  Vj  be 

the  space  spanned  by  •  Set  uj  =  ^  with  coefficients  u^,u;i,u)2,u;3  chosen 

such  that 

f  t  iT * 

/  Wjijjjdxdy  =  |T|  /  —dxdy 
Jr  Jt  a 

I  J,,ujTds  =  \li\  I -^da,  t  =  1,2,3 

•  Jii  Jh  lajj, 

where  [0)/^  is  the  average  of  a  over  the  triauigles  sharing  the  edge  i.e. 

\dTo»x 

ir 


u],  -  /  “  ko«‘ ) 

if/iCdTn] 


(5.3) 

(5.4) 

(5.5) 


The  coefficients  tjo,  •  •  •  satisfying  (5.3)  and  (5.4)  are 


Wi  = 


311^1 


2[a]i4 

20/|r|?^T 

Since  for  any  interior  edge  /  he  values  cjr,  are  the  same  on  both  triangles  T  sharing  the  edge, 
the  function  u  composed  from  u>t  on  every  triangle  is  continuous  and  hence  w  €  Using 

(5.1),  (5.2)  we  get: 

(£*(u,M))^  =  E’Jt 


TeM 

=  E 

rgM 

=  E 

T€M 


<  Cai: 

T^M 


f  UrfjJrdxdy  +  r  E  /  Ji^nds 

Jr  i  (g97- 


where 


C«  =  max 

t^ar.T^M  a  \r 


(5.6) 


Since  w  €  /fo  we  get  from  (3.3) 

(r(u,M))*  <  C.S{e,u;)  <  C.||c||£||u;||£. 

Now  let  for  r  =  {nlf-o 

r'i  -  „,n  It  W^vrf  <i^dy 
^  "iM?.  r?  +  }E;..r? 

where  Vr  €  VV  is  the  function  satisfying  /j-  v^dxdy  =  |T|ro,  Vr«is  =  ll,|r<,  t  =  1, 2, 3.  Therefore 
ijjj  is  the  function  Vt  corresponding  to 

ITinr 

air 


To  = 

Ti  = 


»  =  1,2,3 

W». 
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and  so 


=  0  It  II  Vu,|!*rfxdy 

,  (\n^r  ,  lAl'iMn 

^  ‘^"‘■'nTp-nS-wrJ 

.x/lTm  If /■xlrVl'.fJ?' 


where 

Hence 

and  so,  lining  5.6,  we  obtain 
where 

and  we  get  (4.3). 


<  C^C^Wr, 


r>‘  °'\t 

C.  =  max  -TT-- 
“  icaT.reM  a  / 


\\u\\E<im^Cr)C,£*{u,M) 


riu,M)<Qa{maxCr)\\e\\E 


Qa  =  CaC^  =  niM 


ah 


T€M  a  \t‘ 


6  The  Experimental  Results 

In  this  section  we  will  experimentally  verify  the  reliability  of  veirious  conclusions  that  can  be 
made  from  Theorems  4.1,  4.2  for  the  given  estimator. 

6.1  The  Influence  of  Mesh  Topology 

Theorem  4.1  was  proven  for  general  meshes  without  any  restrictions  on  the  mesh  topology. 
Hence,  it  can  be  expected  that  particular  meshes  with  different  topologies  could  influence  the 
effectivity  index  in  the  same  order  of  magnitude  as  indicated  by  the  theoretical  upper  and 
lower  bounds.  In  addition  the  influence  of  the  angle  has  to  be  visible  (and  could  also  depend 
on  the  topology).  To  check  the  validity  of  this  conclusion  for  a  smooth  solution  we  consider 
the  following  problem: 

For  n  =  (0, 1)  X  (0,1),  let  u  be  the  solution  of 


—  Au  =  /  in  n,  «  =  0  on  T 


(6.1) 


We  choose  /  such  that  u  =  sinxisinxy. 
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mesh  ratio  1/1  mesh  ratio  1/2 

Figure  6.1:  Regula*'  three-direction  mesh 


6.1.1  Three-direction  meshes 

Let  us  consider  the  three-directional  meshes  which  are  obtained  by  uniform  refinement  of  a 
basic  mesh  of  the  type  =  1,2,4,8  as  shown  in  Fig.  6.1. 

In  Table  6.1  we  report  the  effect! vity  index  for  these  meshes.  We  also  report  the  upper  and 
lower  bounds  for  the  effectivity  index  from  Theorems  4.1,  4.2,  when  p  =  2  was  used. 


number  of 

mesh 

elements 

1/1 

1/2 

1/4 

1/8 

8 

0.852 

16 

0.993 

32 

1.070 

1.289 

64 

1.183 

1.804 

128 

1.147 

1.503 

256 

1.252 

2.068 

512 

1.180 

1.571 

1024 

1.277 

2.145 

2048 

1.192 

1.594 

4096 

1.287 

2.169 

upper  bound 

2.34 

2.87 

3.87 

5.38 

lower  bound 

0.29 

0.21 

Table  6.1:  Effectivity  index  for  problem  6.1  (u  =  sinxisiniry)  and  3-direction  meshes 


6.1.2  Criss-cross  meshes 

Let  us  consider  now  the  ‘criss-cross’  meshes  as  shown  in  Fig.  6.2.  The  effectivity  indices  for 
these  meshes  are  given  in  Table  6.2. 


6.1.3  Diamond-shaped  meshes 

Finally,  we  consider  the  ‘diamond-shaped’  meshes  as  shown  in  Fig.  6.3.  The  effectivity  indices 
are  given  in  Table  6.3. 

Comparing  the  tables  we  can  make  the  following  conclusions,  which  are  in  agreement  with 
Theorems  4.1,  4.2: 
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ww 


mwA^^A 


mesh  ratio  1/1 


mesh  ratio  1/2 


Figure  6.2;  Regular  criss-cross  mesh 


number  of 

mesh 

elements 

1/1 

1/2 

1/4 

1/8 

8 

1.350 

16 

0.939 

32 

1.036 

1.049 

64 

1.053 

1.349 

128 

1.077 

1.179 

256 

1.085 

1.518 

512 

1.088 

1.214 

1024 

1.093 

1.562 

2048 

1.090 

1.223 

4096 

1.095 

1.574 

upper  bound 

2.34 

2.87 

3.87 

5.38 

lower  bound 

0.29 

0.21 

0.15 

0.12 

Table  6.3:  EfFectivity  index  for  problem  6.1  (u  =  sinTrxsinjry)  and  ‘criss-cross’  meshes 


mesh  ratio  1/1 


mesh  ratio  1/2 


Figure  6.3:  Regular  ‘diamond-shaped’  mesh 
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1.  For  /t  — »  0  the  effectivity  index  stays  in  the  theoretical  bounds.  In  addition  we  see  that  in 
our  concrete  case  the  effectivity  index  converges  to  a  limiting  value.  This  is  not  surprising, 
because  the  mesh  is  translation  invariant  and  the  solution  can  be  periodically  extended. 
In  this  case  various  results  about  superconvergence  can  be  used  (see  [21]). 

2.  The  effectivity  index  depends  on  the  topology.  The  error  estimator  is  not  asymptotically 
exapt. 

3.  The  dependence  of  the  effectivity  index  on  the  angle  is  clearly  visible. 

4.  The  theoretical  bounds  for  the  effectivity  index  are  relatively  pessimistic.  This  is  because 
these  bounds  are  independent  of  the  topology  and  the  solution.  Better  bounds  in  our 
case  could  be  obtained  by  exploiting  the  superconvergence  effect  (see  [5]). 

6.2  The  Influence  of  the  Solution  Structure 

Theorems  4.1  and  4.2  have  been  shown  for  a  general  solution.  We  expect  that  the  effectivity 
index  may  depend  on  the  solution  structure.  To  check  this  conclusion  consider  the  solution  of 
the  following  problem  on  H  =  (0, 1)  x  (0, 1): 

dxjL 

—  Au  =  /  in  fl,  u  =  0  for  x  =  0  and  x  =  1,  =  0  for  y  =  0  and  y  =  1  (6.2) 

ay 

with  /  such  that  u  =  sinxx.  We  will  use  meshes  as  in  Fig.  6.1  with  the  orientation  shown 
there  (mesh  1/2  and  1/4)  and  with  the  opposite  orientation  (2/1  and  4/1). 


number  of 
elements 

1/1 

1/2 

mesh 

1/4 

2/1 

4/1 

8 

1.013 

16 

1.564 

0.717 

32 

1.107 

2.258 

0.507 

64 

1.598 

0.783 

128 

1.131 

2.273 

0.553 

256 

1.608 

0.800 

512 

1.137 

2.278 

0.566 

1024 

1.611 

0.804 

2048 

1.139 

2.279 

0.569 

4096 

1.612 

0.806 

upper  bound 

2.34 

2.87 

3.87 

2.87 

3.87 

lower  bound 

wm 

Table  6.4:  Effectivity  index  for  problem  6.2  (u  =  sin  xx)  and  3-direction  meshes 

From  Table  6.4  we  see  in  fact  that  the  effectivity  index  depends  on  the  relation  between  the  mesh 
and  the  structure  of  the  solution  (especially  the  relation  between  the  second  derivatives  and  the 
mesh  orientation).  Hence  in  general  we  have  to  deal  with  uncertainties  in  the  effectivity  index. 
(No  correction  factor  could  improve  the  situation.)  We  note  that  our  factor  k  =  0.2684 . . . 
could  be  changed,  but  in  general  (for  all  functions)  the  range  in  the  effectivity  index  could  not 
be  improved.  Again,  we  see  that  Theorems  4.1  and  4.2  reliably  predict  the  outcome  of  the 
numerical  experiments. 
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6.3  The  Influence  of  the  Smoothness  of  the  Solution 


Theorems  4.1  and  4.2  contain  higher  order  terms  which  we  neglected  in  the  reported  theoretical 
bounds  in  Tables  6. 1-6.4.  when  computing  the  bounds.  These  terms  influence  the  convergence 
of  the  effectivity  index  as  A  — »  0.  To  check  the  effect  of  these  terms  we  will  consider  once  more 
problem  6.1,  using  now  the  solution 

*^(^1  y)  =  sin  aitx  sin  axy.  (6.3) 

(We  use  regular  three-direction  meshes,  as  in  Fig.  6.1  with  mesh  ratio  1/1). 


number  of 
elements 

a  =  1 

a  =  2 

tl 

a=  11 

8 

0.852 

0.833 

0.438 

32 

1.070 

0.876 

0.932 

0.379 

128 

1.147 

1.080 

0.887 

0.530 

512 

1.180 

1.161 

1.090 

0.780 

2048 

1.192 

1.187 

1.168 

1.022 

8192 

1.148 

upper  bound 

2.34 

2.34 

2.34 

2.34 

lower  bound 

0.29 

0.29 

0.29 

0.29 

Table  6.5:  Effectivity  index  for  problem  6.3  (u  =  sin  axi  sin  axy)  and  3-direction  meshes 

We  see  that  the  insufficient  smoothness  influences  the  effectivity  index  and  that  the  effec¬ 
tivity  index  can  behave  erratically.  Nevertheless  the  bounds  — where  higher  order  terms  are 
neglected —  are  still  valid.  We  see  in  general  that  the  effectivity  index  decreases  with  decreas¬ 
ing  smoothness  of  the  solution.  The  lower  bound  goes  down  with  increasing  p,  but  the  upper 
bound  is  not  influenced.  Hence  we  can  expect  that  the  effectivity  index  will  go  down  (assuming 
heuristically  that  the  effecitivity  index  moves  as  the  middle  of  the  interval).  The  same  effect 
will  be  seen  later.  Once  more  we  see  good  agreement  with  the  conclusions  stemming  from 
Theorem  4.1  and  that  we  cannot  hope  for  the  existence  of  a  universal  correction  factor. 


6.4  Harmonic  solutions 

The  error  estimator  contains  terms  related  to  the  residual  and  the  jumps,  respectively.  To  isolate 
the  influence  of  the  jump  term,  we  consider  the  following  problem  (again  in  fl  =  (0, 1)  x  (0, 1)): 

• 

-Au  =  0mn,  «  =  0  for  i  =  0,*  =  l  and  y  =  0,  -3- =  y  for  y  =  1  (6.4) 

on 

and  let  u(z,y)  =  sinxzsinhxy.  The  effectivity  index  (for  a  regular  three-direction  mesh  as 
shown  in  Fig.  6.1)  is  given  in  Table  6.6. 

Comparing  Table  6.1  and  6.6  we  see  a  similar  behaviour  influenced  only  mildly  by  the  structure 
of  the  solution.  The  limiting  value  for  the  effectivity  index  seems  to  be  different  for  the  mesh 
1/1.  For  the  other  meshes  this  effect  is  less  visible.  For  these  particular  meshes,  this  might  be 
explained  by  a  superconvergence  phenomenon. 
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number  of 

mesh 

elements 

msM 

mm 

no 

8 

0.841 

16 

0.958 

32 

0.978 

1.269 

64 

1.102 

1.795 

128 

1.040 

1.457 

256 

1.166 

2.042 

512 

1.065 

1.534 

1024 

1.190 

2.136 

2048 

1.071 

1.560 

4096 

1.199 

2.166 

upper  bound 

2.34 

3.87 

5.38 

lower  bound 

0.29 

■owl 

0.15 

■la 

Table  6.6:  EfFectivity  index  for  problem  6.4  (ti(i,y)  =  sinirxsinhjrj/)  and  3-direction  meshes 

6.5  Adaptively  Constructed  Meshes 

So  far,  we  have  restricted  our  experiments  to  uniform  meshes.  We  will  now  consider  nonuniform 
meshes  either  based  on  an  adaptive  procedure  or  randomly  refined.  We  use  problem  6.4  again. 
We  have  constructed  three  sequences  (A  -  C)  of  meshes.  Each  sequence  is  obtained  by  refining 
a  part  of  the  elements  of  the  previous  mesh  according  to  a  refinement  indicator  and  adjusting 
neighbouring  elements.  In  sequences  A  and  B  the  refinement  indicator  was  an  error  estimator 
(with  different  parameters),  in  sequence  C  the  indicator  was  a  random  number.  Typical  meshes 
of  each  sequences  are  shown  in  Fig.  6.4.  The  effectivity  indices  are  given  in  Table  6.7. 


Sequence  A 

Sequence  B 

Sequence  C 

elements  effectivity  index 

elements  effectivity  index 

elements  effectivity  index 

8 

0.841 

8 

0.841 

8 

19 

19 

0.944 

23 

47 

0.970 

38 

0.941 

45 

0.964 

94 

1.008 

67 

96 

1.008 

181 

1.018 

107 

0.949 

200 

0.979 

323 

1.041 

201 

435 

0.990 

593 

1.043 

355 

1.039 

964 

0.972 

1020 

1.060 

607 

1.032 

2129 

1.008 

1773 

1.051 

991 

1.056 

4588 

0.992 

Table  6.7:  Effectivity  indices  for  problem  6.4  and  nonuniform  meshes 


Comparing  Tables  6.7  and  6.6  we  see  that  different  refinement  strategies  lead  to  different 
effectivity  indices,  although  the  difference  is  not  big.  Once  more  we  see  a  r2Lnge  of  the  effectivity 
indices.  We  note  that  for  these  general  meshes  superconvergence  results  cannot  be  applied. 
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Figure  6.4:  Nonuniform  meshes  for  problem  6.4 
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6.6  Effectivity  Indices  for  Problems  with  Unsmooth  Solutions 

So  far,  we  have  considered  problems  with  a  smooth  solution.  Let  us  now  examine  the  behaviour 
of  the  error  estimator  for  solutions  with  singularities.  Theorems  4.1  and  4.2  indicate  a  decrease 
of  the  effectivity  index,  when  the  singularity  of  the  exact  solution  grows,  and  also  that  this 
effect  diminishes  when  adaptive  meshes  are  used. 

To  this  purpose  we  consider  the  domain  fl  =  (—1,1)  x  (—1,1),  partitioned  into  4  square 
subdomains  as  shown  in  Fig.  6.5.  Let  u  be  the  solution  of 

dti 

-  V(a(i,j/)Vu)  =  0  in  n,  a(i,y)^  =  p(i,  j/)  on  F  (6.5) 

where  a(x,y)  is  constant  in  each  quadrant  (a(x,y)  =  1  or  a(x,y)  =  a,  see  Fig.  6.5).  a  is  chosen 
so  that 

u(x,  y)  =  r"(cj  cos(a;9)  +  s,  sin(Qt¥>))  (6.6) 

where  are  polar  coordinates  with  center  in  the  origin.  For  a  =  0.5  and  a  =  0.1  we  get 
a  =  3  +  «  5.828 . . .  and  a  ss  166.447 . . .,  respectively. 


2 

a(x,y)  =  a 

1 

a(x,y)=  1 

3 

a(x,y)=  1 

4 

o(x, y)  =  0 

Figure  6.5:  Coefficient  a  for  problem  6.5 


Table  6.8  shows  the  effectivity  index  for  uniform  meshes.  (In  each  quadrant  the  meshes  are 
uniform  three-direction  meshes  of  mesh  ratio  1/1  as  in  Fig.  6.1.  The  orientation  of  the  meshes 
is  changed  in  subdomuns  2  and  4  to  preserve  the  symmetry  of  the  solution.) 


number  of 

elements 

a  =  0.5  a  =  0.1 

32 

128 

512 

0.719  0.536 

8192 

0.722 

Table  6.8:  Effectivity  index  for  problem  6.5  (u  =  r‘*f{ip))  and  3-direction  meshes 

Although  the  estimator  in  Theorem  4.1  includes  the  ratio  of  a,  it  influences  only  elements  on 
the  interfaces  and  hence  this  factor  can  be  neglected.  The  decrease  of  the  effectivity  index  is 
clearly  related  to  the  necessity  to  use  a  higher  degree  p  in  4.2  for  the  estimate.  Hence  the 
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unsmoothness  of  the  solution  decreases  the  effectivity  index.  This  has  also  been  observed  in 
Table  6.5. 


We  have  seen  in  Table  6.7,  that  nonuniformity  of  the  mesh  does  not  influence  the  effectivity 
index  too  much.  On  the  other  hand  an  adaptive  procedure  refines  the  mesh  in  the  places  where 
the  solution  in  unsmooth  and  hence  a  lower  p  can  be  used.  We  have  created  two  sequences  of 
adaptively  refined  meshes  (A  and  B).  The  meshes  B1  and  B2  are  less  refined  than  the  meshes 
A1  and  A2.  The  refinement  in  meshes  B  is  more  concentrated  around  the  origin  than  in  meshes 
A.  Examples  of  each  sequences  are  displayed  in  Fig.  6.6.  If  the  mesh  is  properly  refined  ,  a 
lower  p  can  be  used  and  the  effectivity  index  improves.  Table  6.9  confirms  this  conclusion  (for 
a  =  0.5). 
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Figure  6.6:  Nonuniform  meshes  for  problem  6.5 


The  effect  of  the  unsmoothness  of  the  solution  can  also  be  observed  if  the  singular  part  of  the 
solution  is  removed.  So  we  consider  problem  6.5  in  the  domain  fi  =  —  [—0.5, 0.5]  x  [—0.5, 0.5] 

and  prescribe  =  <;  on  dCl.  Table  6.10  shows  the  effectivity  indices  for  this  problem  using 
uniform  meshes  as  for  Table  6.8.  Again,  a  =  0.5. 

We  note  that  the  effectivity  index  is  now  very  similar  to  problems  6.1,  6.2. 
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Sequence  A 

Sequence  B 

elements  effectivity  index 

elements  effectivity  index 

32 

0.658 

32 

87 

64 

194 

0.738 

96 

0.787 

370 

168 

0.830 

725 

0.787 

244 

0.880 

1297 

0.807 

363 

0.918 

2346 

0.828 

485 

0.960 

604 

0.959 

768 

0.975 

965 

0.995 

Table  6.9:  EfFectivity  indices  for  problem  6.5  (u  =  and  adapted  meshes 


number  of 
elements 

a  =  0.5 

24 

0.964 

96 

1.027 

384 

1.059 

1536 

1.070 

6144 

1.073 

Table  6.10:  Effectivity  index  for  problem  6.5  (u  =  on 


n)  wd  3-direction  meshes 
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6.7  The  Oscillation  Effect 


It  is  often  assumed,  that  the  error  of  the  interpolant  is  close  to  the  error  of  the  finite  element 
method.  To  demonstrate  that  this  assumption  is  in  gener<d  incorrect,  we  show  that  the  error 
of  the  finite  element  method  could  oscillate.  To  this  end  we  show  in  Fig.  6.7  the  typical 
error  behaviour  for  the  problem  discussed  in  section  6.6  with  uniform  meshes  and  a  =  0.5. 
Displayed  are  two  patches  of  regular  three-direction  meshes  at  the  same  location — far  from  the 
singularity — in  D.  Fig.  6.7  also  shows  the  error  indicators  (in  parentheses).  In  Fig.  6.8  we 
illustrate  graphically  the  error  behaviour  for  the  entire  mesh  where  the  mesh  has  363  elements. 
(An  element  is  black  if  the  energy  norm  error  in  this  element  exceeds  a  certain  threshold.  Here, 
the  threshold  is  4.0  x  10~^.) 
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Figure  6.7:  Oscillation  of  the  finite  element  error 

We  see  that  the  error  indicator  shows  approximately  the  average  error  for  adjacent  elements 
and  hence  does  not  exhibit  the  oscillation  of  the  error  of  the  finite  element  solution.  Note  that 
the  interpolamt  does  not  show  the  oscillation.  This  oscillation  effect  is  caused  by  the  pollution 
of  the  error  stemming  from  the  singularity  of  strength  If  we  enforce  a  smooth  solution  by 
removing  the  center  part  of  the  domain,  as  in  problem  6.5,  the  results  given  in  Fig.  6.9  show 
that  the  oscillation  disappears.  (We  use  elements  from  the  same  location  as  in  Fig.  6.7;  the 
numbers  are  therefore  directly  comparable.) 

The  oscillation  effect  is  weaker,  if  a  ^  0.5  or  the  mesh  is  constructed  adaptively.  We  show  in 
Fig.  6.10  the  local  results  for  a  =  .8  and  in  Fig.  6.11  the  results  for  o  =  .5  and  adaptively 
refined  meshes.  In  Fig.  6.12  we  show  the  error  behaviour  for  an  adaptively  refined  mesh  and 
the  threshold  2.0  x  10“*.  In  Section  7  we  will  examine  this  effect. 


7  Analysis  of  the  Oscillation  Phenomenon 

The  oscillation  phenomenon  we  have  observed  in  Section  6.7  occurs  also  in  the  case  of  the 
following  problem  on  fl  =  [—1, 1]  x  [0, 1) 

-Au  =  0  in  n,  u  =  0  on  To  =  {(i.y)  i  0  <  x  <  l,y  =  0),  ~  ^  r;v  =  r-rD 
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small  element  side:  1/8  small  element  side:  1/16 

values  in  (..):  error  estimator 

Figure  6.9:  Finite  element  error  and  estimators  for  a  problem  without  singularity 
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Figure  6.10:  Finite  element  error  and  estimators  for  problem  with  singularity  (a  =  .8,  uniform 
mesh) 
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Figure  6.11:  Finite  element  error  and  estimators  for  problem  with  singularity  (a  =  .5,  adap* 
tively  refined  mesh) 
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Figure  6.12:  Finite  element  error  for  adapted  mesh  (a  =  0.5) 
with  g  selected  so  that 

u  =  r*/^sin(ivj)  =  Rez*^* 

with  z  =  X  +  iy. 

We  consider  now  a  uniform  3-direction  mesh  of  the  type  shown  in  Fig.  6.1.  Then  the  finite 
element  solution  has  essentially  the  form 

«s(x,y)  =  u^(i,y)-hG(x,y)-t-i?(i,y),  (7.1) 

where  .)||/f‘  ^  o(^)  G(i,y)  is  a  piecewise  line.ar  function,  which  coincides  in  the  nodal 
points  of  the  mesh  with  the  function 

G*(i,y)  =  sin  =  y4APez“’'*.  (7.2) 

The  function  G’{x,y)  is  the  adjoint  singularity  function.  It  describes  the  pollution  error,  which 
is  of  the  same  order  (in  ||.||e)  as  the  error  of  the  interpolant.  The  error — measured  in  the 
energy  norm — is  locally  (i.e.  on  single  elements)  of  order  and  is  governed  by  the  second 
derivatives  of  u.  The  contributions  of  G  to  the  energy  err  r  are  of  the  same  order,  but  are 
governed  by  the  first  derivative  of  G.  Now  observe  that  in  our  case  both  errors  are  governed  by 
the  second  derivative  of  Rez^^^.  This  leads  to  a  cancellation  and  consequently  to  the  oscillation 
phenomenon. 

To  demonstrate  this,  we  consider  a  patch  of  two  triangles  as  in  Fig.  7.1  with  local  coordinates 
f ,  Tf  with  the  origin  in  the  lower  left  vertex  of  the  patch. 

Using  (7.1),  (7.2)  we  can  assume  that  on  the  patch: 

u  =  a{(^  —  rj^)  +  b^Tj  -J-  linear  functions  +  h.o.t. 

G  =  (aa^  -I-  abri)h  constant  term  +  h.o.t. 
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Figure  7.1:  Mesh  for  local  analysis  of  errors 
where  a  is  related  to  A  in  (7.2).  Neglecting  R  in  7.1,  we  get  by  simple  computations 


Ikllli  _  +  2a  +  3q^)  +  g" 

INIIs  (a2  +  62)(l-2o  +  3a2)  +  a* 

(||e||ei  is  the  energy  norm  of  the  error  in  element  t)  and  we  see  that  we  have  in  the  case  without 
pollution  (a  =:  0) 

Mli-1 

llelll, 

while  the  ratio  is  ^  1  in  the  presence  of  the  pollution  error.  This  explains  our  observation  in 
Section  6.7. 


Let  us  note  that  the  said  effect  occurs  when  the  pollution  error  and  the  error  of  the  inter- 
polant  have  the  same  strength.  This  is  in  agreement  with  our  results  in  Section  6.7,  where  the 
oscillation  phenomenon  was  only  visible  for  the  singularity  a  =  .5. 


8  Conclusions 


Summarizing  the  results  that  we  have  presented  in  the  previous  sections  we  draw  the  following 
conclusions. 

1.  The  estimator  performs  as  expected  from  its  properties  given  in  Theorems  4.1  and  4.2. 

2.  The  effectivity  index  depends  on  the  mesh,  especially  on  the  angle  of  the  triangles.  Hence 
meshes  without  small  angles  are  preferable.  The  notion  of  the  angle  depends  on  the 
differential  operator. 

3.  The  effectivity  index  is  not  too  sensitive  to  the  topology  of  the  mesh  (except  for  the 
minimal  angle). 

4.  The  effectivity  index  can  be  larger  or  smaller  than  1  depending  on  the  character  of  the 
solution  as  predicted  in  Theorems  4.1  and  4.2. 

5.  If  the  solution  has  singular  behaviour,  the  quality  of  the  error  estimator  deteriorates.  The 
deterioration  can  be  avoided  by  adaptive  refinement,  as  follows  from  Theorem  4.1. 
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